function npde(tau);
close all

kmax=500;
n=50;
theta=0;
a=1;
b=1;
h=1/(n+1)
q=tau/h^2
M0=1;

r=zeros(1,n);
r(1,1)=-2;
r(1,2)=1;

Q=toeplitz(r);

s=zeros(1,n);
s(1,1)=0;
s(1,2)=-1;
r(1,1)=0;
r(1,2)=1;

D=toeplitz(s,r);
I=eye(n,n);

%initial vector

f=inline('sin(x*pi)+3*sin(2*x*pi)-10*sin(5*x*pi)');

v=zeros(n,1);
for i=1:n
   v(i)=f(i*h);   
end;

   w(1)=0; w(n+2)=0;
   w(2:n+1)=v;
   
   plot([0:h:1],w)   
   
for k=1:kmax
      
   for i=1:n
      fdiag(i)=sin(v(i))*v(i);
   end;
   
   F=diag(fdiag);
   
   A1=I-theta*((a*tau/h^2)*Q+(b*tau/2/h)*D+tau*F); A1=sparse(A1);
   A2=I+(1-theta)*((a*tau/h^2)*Q+(b*tau/2/h)*D+tau*F); A1=sparse(A1);
   
   v=A1\A2*v;
   
   w(1)=0; w(n+2)=0;
   w(2:n+1)=v;
   
   
   plot([0:h:1],w)
   
      pause(0.1)
end;

