%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% heat flow solved by use of Heun - RK
%
% VHS Nov. 2003
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf;

% h = time-grid size, n = number of time steps

%h = 0.01; n=25;
h = 0.00495; n = 51;

[X,Y] = meshgrid(1:1:n,0:0.1:1);

% initial data
for j = 1:1:11
    Z(j,1) = sin((j-1)*pi/10);
end

% method of Heun
for k = 2:1:n
    Z(1,k) = 0;
    Z(11,k) = 0;
    for i = 2:1:10
        f0(i) = 100*(Z(i+1,k-1)-2*Z(i,k-1)+Z(i-1,k-1));
        u1(i) = Z(i,k-1) + h*f0(i);
    end
    u1(1) = 0; u1(11)=0;
    for i = 2:1:10
        f1(i) = 100*(u1(i+1)-2*u1(i)+u1(i-1));
    end
    for i = 2:1:10
        Z(i,k) = Z(i,k-1) + h*(0.5*f1(i)+0.5*f0(i));
    end
end

% 3D-plot
meshc(X,Y,Z);