Applications - A Guide to MATLAB for Beginners and Experienced Users - page 190

where s j = k j t / ( x ) 2 . In the following M-file, which we called heatvc.m ,
we modify our previous M-file to incorporate this iteration.
function u = heatvc(k, x, t, init, bdry)
% Solve the 1D heat equation with variable coefficient k on
% the rectangle described by vectors x and t with
% u(x, t(1)) = init and Dirichlet boundary conditions
% u(x(1), t) = bdry(1), u(x(end), t) = bdry(2).
J = length(x);
N = length(t);
dx = mean(diff(x));
dt = mean(diff(t));
s = k*dt/dxˆ2;
u = zeros(N,J);
u(1,:) = init;
for n = 1:N-1
u(n+1, 2:J-1) = s(2:J-1).*(u(n, 3:J) + u(n, 1:J-2)) + ...
(1 - 2*s(2:J-1)).*u(n,2:J-1) + ...
0.25*(s(3:J) - s(1:J-2)).*(u(n, 3:J) - u(n, 1:J-2));
u(n+1, 1) = bdry(1);
u(n+1, J) = bdry(2);
end
Notice that k is now assumed to be a vector withthe same lengthas x and
that as a result so is s . This in turn requires that we use vectorized
multiplication in the main iteration, which we have now split into three
lines.
Let's use this M-file to solve the one-dimensional variable-coefficient heat
equation withthe same boundary and initial conditions as before, using
k ( x ) = 1 + ( x / 5) 2 . Since the maximum value of k is 2, we can use the same
values of t and x as before.
kvals=1+(xvals/5).ˆ2;
uvals = heatvc(kvals, xvals, tvals, init, [15 25]);
surf(xvals, tvals, uvals)
xlabel x; ylabel t; zlabel u