I think I can change the for loops with seq but I don't know how to do this for the boundary and initial conditions.

Example: I think I can replace

for i from 0 to nxpoints do
u[i,0]:=sin(Pi*i*h)
end do:

With something along the lines of

seq(u[i,0]:=sin(Pi*i*h), i=0..nxpoints)
end do;

Although when I do this I get errors further up the code.

Thank you for your reply. I have replaced

Initial condition
for i from 0 to nxpoints do
u[i,0]:=sin(Pi*i*h)
end do:
solution:=[[0.,seq(u[i,0],i=0..nxpoints)]]:
printf(format,op(solution[1])):

with

Initial condition
seq(assign( u[i,0]=sin(Pi*i*h) ),i=0..nxpoints);
end do:
solution:=[[0.,seq(u[i,0],i=0..nxpoints)]]:
printf(format,op(solution[1])):

After the first plot command it gives me and error and pulls me back to the seq. I can't seen a reason for the error, the code up until the error is below. Any help again is greatly appreciated.

restart;
h:=0.1: k:=0.001: r:=k/h^2:
xl:=0: xr:=1: # left and right ends of interval in x
nxpoints:=round((xr-xl)/h):
ntsteps:=5: # total number of time steps
iprint:=1: # print output every iprint time steps
printf(cat(" x,t ","%7.3f"$(nxpoints+1),"\n"),seq(i*h,i=0..nxpoints)):
format:=cat("%7.3f"$(nxpoints+2),"\n"): # format for printing output
Boundary conditions
for j to ntsteps+1 do
u[0,j]:=0: u[nxpoints,j]:=0
end do:
Initial condition
seq(assign(u[i,0]=sin(Pi*i*h)),i=0..nxpoints);
end do:
solution:=[[0.,seq(u[i,0],i=0..nxpoints)]]:
printf(format,op(solution[1])):
Time steps
for j from 0 to ntsteps-1 do
for i from 1 to nxpoints-1 do
u[i,j+1]:=r*u[i-1,j]+(1-2*r)*u[i,j]+r*u[i+1,j]
end do;
solution:=[op(solution),[(j+1)*k,seq(u[i,j+1],i=0..nxpoints)]]:
if floor((j+1)/iprint)=(j+1)/iprint then printf(format,op(solution[j+2])) end if;
end do:
Plots
for j from 1 by iprint to ntsteps+1 do
plot||j:=plot([seq([i*h,solution[j,i+2]],i=0..nxpoints)],labels=[x,u]):
end do:
plots[display]({seq(plot||(iprint*j+1),j=0..ntsteps/iprint)});
x,t 0.000 0.100 0.200 0.300 0.400 0.500 0.600 0.700 0.800 0.900 1.000
Error, reserved word `end` unexpected

This is my code in its current state.

restart;
h:=0.1: k:=0.001: r:=k/h^2:
xl:=0: xr:=1: # left and right ends of interval in x
nxpoints:=round((1/2)*(xr-xl)/h):
ntsteps:=5: # total number of time steps
iprint:=1: # print output every iprint time steps
printf(cat(" x,t ","%7.3f"$(nxpoints+1),"\n"),seq(i*h,i=0..nxpoints)):
format:=cat("%7.3f"$(nxpoints+2),"\n"): # format for printing output
Boundary conditions
for j to ntsteps+1 do
u[0,j]:=0: u[nxpoints,j]:=0
end do:
Initial condition
seq(assign(u[i,0]=sin(i*h)),i=0..nxpoints);
solution:=[[0.,seq(u[i,0],i=0..nxpoints)]]:
printf(format,op(solution[1])):
Time steps
for j from 0 to ntsteps-1 do
for i from 1 to nxpoints-1 do
u[i,j+1]:=r*u[i-1,j]+(1-2*r)*u[i,j]+r*u[i+1,j]
end do;
i=nxpoints: u[i,j+1]:=2*r*u[i-1,j]+(1-2*r)*u[i,j]; # at the point x=0.500, taking account of symmetry
for i from nxpoints to 1 do
u[i,j]=u[1-i,j]
end do:
solution:=[op(solution),[(j+1)*k,seq(u[i,j+1],i=0..nxpoints)]]:
if (j+1) mod iprint = 0 then printf(format,op(solution[j+2])) end if;
end do:
Plots
for j from 1 to nops(solution) do
plot||j:=plot([seq([i*h,solution[j,i+2]],i=0..nxpoints)],labels=[x,u]):
end do:
plots[display]({seq(plot||(iprint*j+1),j=0..ntsteps/iprint)});
x,t 0.000 0.100 0.200 0.300 0.400 0.500
0.000 0.000 0.100 0.199 0.296 0.389 0.479
0.001 0.000 0.100 0.198 0.295 0.389 0.461
0.002 0.000 0.100 0.198 0.295 0.387 0.447
0.003 0.000 0.100 0.198 0.294 0.384 0.435
0.004 0.000 0.099 0.198 0.294 0.380 0.425
0.005 0.000 0.099 0.198 0.293 0.376 0.416

Comparison with exact solution
plottime:=0.005:
uexact:=(x,t)->4*cos(1/2)*sum(((-1)^n)/((2*n+1)^(2)*Pi^(2)-1)*exp(-(2*n+1)^(2)*Pi^(2)*t)*sin((2*n+1)*Pi*x),n=0..100):
pexact:=plot(uexact(x,plottime),x=0..xr,labels=[x,u],linestyle=DASH):
plots[display]({pexact,plot||(floor(plottime/k)+1)});

I have changed my definition of nxpoints at the top in a crude attempt to take advantage of the symmetry of the problem. This has cut down on the computations needed but I have not been able to plot the whole graph so this symettry isn't being carried through. I have tried to enter the extra values needed using this code:

for i from nxpoints to 1 do
u[i,j]=u[1-i,j]
end do:

It seems rather ineffective, it doesn't change any of the output or graphs and also doesn't result in any errors. Plotting the whole graph is the last piece of the puzzle, any changes to code after getting this sorted would be more for efficency then actually getting the answer.

Doug, I used your suggestion for X mod N = 0 and it has worked nicely. I tried to use nested seq's but failed miserably.

Finally, I tried to replace

for j from 0 to ntsteps-1 do
for i from 1 to nxpoints-1 do
u[i,j+1]:=r*u[i-1,j]+(1-2*r)*u[i,j]+r*u[i+1,j]
end do;

with

seq(assign(u[i,j+1]=r*u[i-1,j]+(1-2*r)*u[i,j]+r*u[i+1,j]),j=0..nxsteps-1,i=1..nxpoints-1);

it didn't work and threw up many errors.

I have put some code in to try and give results for x between 0.6 and 1 (the 6th to 10th steps), I have added

for i from nxpoints+1 to 2*nxpoints do
u[i,j]=u[1-i,j]
end do;

Into my timestep loop immediately before

i=nxpoints: u[i,j+1]:=2*r*u[i-1,j]+(1-2*r)*u[i,j];

The problem is my results are only printing up to x=0.5 leaving half the solution out. This problem follows into my graphs and I only get a comparison up to x=0.5. I think I need to extend the output up to the correct number of points but I'm finding this difficult due to my definition of nxpoints at the top of the code.

Any ideas?

Doug, thanks for your help. I have finally got a solution that works.