> restart;
> # HEAT EQUATION BACKWARD-DIFFERENCE ALGORITHM 12.2
> #
> # To approximate the solution to the parabolic partial-differential
> # equation subject to the boundary conditions
> #                u(0,t) = u(l,t) = 0, 0 < t < T = max t,
> # and the initial conditions
> #                u(x,0) = F(x), 0 <= x <= l:
> #
> # INPUT:   endpoint l; maximum time T; constant ALPHA; integers m, N.
> #
> # OUTPUT:  approximations W(I,J) to u(x(I),t(J)) for each
> #          I = 1, ..., m-1 and J = 1, ..., N.
> alg122 := proc() local F, OK, FX, FT, ALPHA, M, N, M1, M2, N1, H, K, VV, I, W, L, U, J, T, Z, I1, FLAG, NAME, OUP, X;
> printf(`This is the Backward-Difference Method for Heat Equation.\n`);
> printf(`Input the function F(X) in terms of x.\n`);
> printf(`For example:  sin(3.141592654*x)\n`);
> F := scanf(`%a`)[1];
> F := unapply(F,x);
> printf(`The lefthand endpoint on the X-axis is 0.\n`);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the righthand endpoint on the X-axis.\n`);
> FX := scanf(`%f`)[1];
> if FX <= 0 then
> printf(`Must be positive number.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the maximum value of the time variable T.\n`);
> FT := scanf(`%f`)[1];
> if FT <= 0 then
> printf(`Must be positive number.\n`);
> else
> OK := TRUE;
> fi;
> od;
> printf(`Input the constant alpha.\n`);
> ALPHA := scanf(`%f`)[1];
> OK := FALSE;
> while OK = FALSE do
> printf(`Input integer m = number of intervals on X-axis\n`);
> printf(`and N = number of time intervals - separated by a blank.\n`);
> printf(`Note that m must be 3 or larger.\n`);
> M := scanf(`%d`)[1];
> N := scanf(`%d`)[1];
> if M <= 2 or  N <= 0 then
> printf(`Numbers are not within correct range.\n`);
> else
> OK := TRUE;
> fi;
> od;
> if OK = TRUE then
> M1 := M-1;
> M2 := M-2;
> N1 := N-1;
> # Step 1
> H := FX/M;
> K := FT/N;
> VV := ALPHA*ALPHA*K/(H*H);
> # Step 2
> for I from 1 to M1 do
> W[I-1] := F(I*H);
> od;
> # Step 3
> # Steps 3 - 11 solve a tridiagonal linear system using Algorithm 6.7
> L[0] := 1+2*VV;
> U[0] := -VV/L[0];
> # Step 4
> for I from 2 to M2 do
> L[I-1] := 1+2*VV+VV*U[I-2];
> U[I-1] := -VV/L[I-1];
> od;
> # Step 5
> L[M1-1] := 1+2*VV+VV*U[M2-1];
> # Step 6
> for J from 1 to N do
> # Step 7
> # Current t
> T := J*K;
> Z[0] := W[0]/L[0];
> # Step 8
> for I from 2 to M1 do
> Z[I-1] := (W[I-1]+VV*Z[I-2])/L[I-1];
> od;
> # Step 9
> W[M1-1] := evalf(Z[M1-1]);
> # Step 10
> for I1 from 1 to M2 do
> I := M2-I1+1;
> W[I-1] := evalf(Z[I-1]-U[I-1]*W[I]);
> od;
> od;
> # Step 11
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:  A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> fprintf(OUP, `BACKWARD-DIFFERENCE METHOD\n\n`);
> fprintf(OUP, `  I         X(I)     W(X(I),%12.6e)\n`, FT);
> for I from 1 to M1 do
> X := I*H;
> fprintf(OUP, `%3d %11.8f %14.8f\n`, I, X, W[I-1]);
> od;
> fi;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> RETURN(0);
> end;
> alg122();