////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//solution to computational problem of hw#6 using backward differencing in 
//space and forward euler in time
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

xinitial=input("Enter x initial>");
xfinal=input("Enter x final>");
tinitial=input("Enter t initial>");
tfinal=input("Enter t final>");

J=input("Enter number of space steps>");
dt=input("Enter time step>");

dx=(xfinal-xinitial)/J;
N=(tfinal-tinitial)/dt;

lambda=dt/dx;

uinitial=zeros(J+1,1);
uk=zeros(J+1,1);
ukt1=zeros(J+1,1);
ukt1=zeros(J+1,1);
ukt2=zeros(J+1,1);
ukt3=zeros(J+1,1);
x=linspace(xinitial,xfinal,J+1)';
t=zeros(tinitial,tfinal,N+1);

uinitial(:)=init(x(:));
uk(:)=uinitial(:);

for n=1:N
	ukt1=dt*Fsystem(uk,t(n));
	ukt2=dt*Fsystem(uk+(1/2)*ukt1,t(n)+(1/2)*dt);

	ukp1=uk+ukt2;
	uk=ukp1;
end

plot2d(x, [uinitial,uk], leg="exact solution@numerical solution");