n = 100;
m = 1000;
x = linspace(0,1,n);
h = x(2)-x(1);
dt = 0.1*h^2;
u(1:n,1) = sin(pi*x);
u(1,2:m) = 0;
u(n,2:m) = 0;
for j=1:m-1
for i=2:n-1
u(i,j+1) = u(i,j) + dt/h^2*(u(i+1,j)-2*u(i,j)+u(i-1,j));
end
end
u
% 열 방정식에서 u(0,t)=0, u(1,t)=0, u(x,0)=sin(pi*x) %
% 위의 boundary condition를 가질때 solution이다%
% 실제로 코드 내에서 u(a,b)의 값은 u(a/99,(b-1)*dt)의 값이다%