ex36ch2.m

these codes are for solving OED with matlab

function ex36ch2
% Quench front problem
global J h u_c A N
A = 2.0d5;
u_c = 0.5;
N = 50;
% Spatial grid:
h = 1 / (N + 1);
x = h * (1:N);
% Initial solution:
u0 = zeros(N,1);
for i = 1:N
    if x(i) >= 0.25
        u0(i) = 1;
    else
        if x(i) > 0.1
            u0(i) = (x(i) - 0.1) / 0.15;
        end
    end
end
% Jacobian matrix:
e = ones(N,1);
J = spdiags([e -2*e e], -1:1, N, N) / h^2;
J(N,N-1) = 2 / h^2;
% ODE solution:
options = odeset('Jacobian', J, 'Stats', 'on');
tic
% Nonstiff solution:
[t,u] =  ode23(@f, [0:0.001:0.006], u0, options);
% Stiff solution:
% [t,u] = ode15s(@f, [0:0.001:0.006], u0, options);
toc
% Add the boundary conditions (extra column on each end):
x = [0 x 1];
npts = length(t);
u = [zeros(npts,1) u zeros(npts,1)];
% Right boundary condition:
u(:,N+2) = u(:,N+1);
% Plot the solution:
surf(x,t,u)
%===================================================================
function dudt = f(t,u)
global J h u_c A N
% g = zeros(N,1);
g = -A * (u .* any(u<=u_c,2));
dudt = J * u + g;