Doug, I have followed all your advice and have completed my code except for one final thing. The original problem was U[T]=2U[XX] with X between 0 and 4, T greater than 0. Subject to the initial condition U(X,0)=0 with X greater than or equal to 0 and less than or equal to 4 and the boundary conditions U[X](0,T)=1 and U[X](4,T)=2, T greater than 0. Calculate up to T=2. I have had to rescale this problem to the standardised form by substituting U=u, X=4x and T=8t. My output and both graphs are in terms of x,t but need to be in terms of X,T. I have played around with some of the labels but cannot get my output and therefor my graphs in the correct terms. Working with the standardised form I have produced the following code: 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:=250: # total number of time steps iprint:=25: # 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 Initial condition for i from 0 to nxpoints do u[i,0]:=0 end do: solution:=[[0.,seq(u[i,0],i=0..nxpoints)]]: printf(format,op(solution[1])): Time steps for j from 0 to ntsteps do u[0,j+1]:=(1-(2*r))*u[0,j] + u[1,j]*2*r - 8*h*r: for i from 1 to nxpoints-1 do u[i,j+1]:=u[i,j] + r*(u[i-1,j] - 2*u[i,j] + u[i+1,j]): end do; u[nxpoints,j+1]:= u[nxpoints,j]*(1-2*r) + 2*r*u[nxpoints-1,j] + 16*h*r: 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)}); Comparison with exact solution plottime:=0.025: uexact:=(X,T)->1/8*(4*T+X^2)+X-8/3-8/Pi^2*sum((2*(-1)^n-1)*exp(-n^2*Pi^2*T/8)*cos(1/4*n*Pi*X)/n^2,n=1..100); pexact:=plot(uexact(x,plottime),x=0..xr,labels=[X,u],linestyle=DASH, color=blue): plots[display]({pexact,plot||(floor(plottime/k)+1)}); Any suggestions are greatly appreciated.

This is the code I have so far for the initial conditions and boundary conditions. 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 Initial Condition for i from 0 to nxpoints do u[i,0]:=sin(i*h) end do: Boundary Condition for j from 1 to ntsteps+1 do u[0,j]:=4: u[nxpoints,j]:=8*h+u[nxpoints-1,j] end do: solution:=[[0.,seq(u[i,0],i=0..nxpoints)]]: printf(format,op(solution[1])): I'm pretty sure the boundary conditions are wrong. Sorry for making you dumb this down even more, this is all completely new to me.

I have made a mistake in the previous post. The boundary conditions are u[x](0,t)=4 and u[x](1,t)=8 I understand the bit about C but I don't understand where (u[N,j]-u[N-1,j])/h=4 comes from. Is N=i and can I put (u[N,j]-u[N-1,j])/h=4 in as a sequence? Something along the lines of seq ( assign ( ( u[N,j] - u[N-1,j] ) / h ) = 4), N=0..nxpoints ); Why do we change N and not j. We know at the boundary that N=0 and N=1.

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

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?

Sorry, I spent 10 minutes putting the above matrices into a more readable format and it doesn't seem to have worked. The general form of the first matrix is a[i]x[i-1] + b[i]x[i] + c[i]x[i+1] = d[i] and the general form for the 2nd matrix is x[i] = g[i]x[i+1] = h[i] Where g[1] = c[1]/b[1], g[i] = c[i]/(b[i] − a[i]g[i−1]), i = 2, 3, . . . , (N − 1) h[1] = d[1]/b[1], h[i] = (d[i] − a[i]h[i−1])/(b[i] − a[i]g[i−1]), i = 2, 3, . . . ,N

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.

Oops, cut that a bit short. Also, with regards to the arrays, we have been banned from using the maple array. Not sure of the reasoning behind it. Thanks again.

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

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.