function ch2ex8 % Galerkin scheme as in Fletcher, p. 112. global J h N = 20; h = 1/(N+1); x = h*(1:N); v0 = x + sin(pi*x); e = ones(N,1); M = spdiags([e 4*e e],-1:1,N,N)/6; J = spdiags([e -2*e e],-1:1,N,N)/h^2; options = odeset('Mass',M,'Jacobian',J); [t,v] = ode15s(@f,[0 0.5],v0,options); % Add the boundary values: x = [0 x 1]; npts = length(t); v = [zeros(npts,1) v ones(npts,1)]; colormap(gray) surf(x,t,v) % generate ch2fig7.eps from o view of this. %========================================= function dvdt = f(t,v) global J h dvdt = J*v; dvdt(end) = dvdt(end) + 1/h^2;