clc
clear
M=80;%x
N=80;%T
L=1;T=1;
h=L/M;%空间步长
tao=T/N;%时间步长
x=0:h:L;
t=0:tao:T;
r=L*(tao/h)^2
f=@(x,t) 0;
%% 构造边值矩阵 
u=zeros(N+1,M+1);
u(1,:)=sin(4*pi*x); 
u(:,1)=0;
u(:,M+1)=0;
%% 方法一 偏心差分 求出网格中的第二行
% for j=2:M
%     u(2,j)=u(1,j)+tao*sin(8*pi*x(j));
% end
%% 方法二 中心差分 求出网格中的第二行
for j=2:M
    u(2,j)=r*(u(1,j-1)+u(1,j+1))/2+(1-r)*u(1,j)+tao*sin(8*pi*x(j));
end
%% 数值解
for i=2:N
    j=2;
    while j<=M;
        u(i+1,j)=r*(u(i,j+1)+u(i,j-1))+2*(1-r^2)*u(i,j)-u(i-1,j)+tao^2*f(x(j),t(i));
        j=j+1;
    end
end
%% 解析解
for j=1:M+1
    for i=1:N+1
    Accurate(i,j)=cos(4*pi*t(i))*sin(4*pi*x(j))+(sin(8*pi*t(i))*sin(8*pi*x(j)))/(8*pi);
    end
end
%% 误差
err=max(max((abs(Accurate-u))))
% %% 画图
% [xx,tt]=meshgrid(x,t);
% figure 
% surf(xx,tt,u)
% xlabel('x');ylabel('t');zlabel('u');title('数值解')

% figure
% surf(xx,tt,Accurate)
% xlabel('x');ylabel('t');zlabel('Accurate');title('解析解')