Applications  
Part I

Sparse equations solving

Mathematica provides the function Solve to solve systems of equations.   When the system is linear in the variables, it can be represented and solved in a matrix form; this functionality is available from LinearSolve. 
 
This uses Solve to solve this system.

Solve[ {3. x  + 2. y == 5., 5. x + 3. y == 1.}, {x, y}]

This uses LinearSolve to solve the same system represented as a matrix problem.

LinearSolve[ {{3., 2.}, {5., 3.}}, {5., 1.}]

Under appropriate conditions, Solve calls the same LU  decomposition routines that LinearSolve uses. These conditions are that  the set of equations is linear, and that the coefficients consist only of machine numbers.  This technique for solving linear  systems is suited to dense matrices where the number of zero entries is small.  For these dense systems, the representation of the  equations is more compact for LinearSolve than for Solve.   
 
Alternatively if the system is sparse,  i.e., if it has a large number of zero entries, the Solve representation becomes much more efficient.  Not only is there an efficient representation for these sparse systems, but there are specialized algorithms to solve them.  Such algorithms have been added to Version 2.2, to be called automatically from Solve when the coefficients are real or complex machine numbers and the system is sparse. 
 
This sets up the equations.

(n = 4000; 
 vars = Table[x[i], {i, n}]; 
 eqns = Join[{-2. x[1] + x[2] == 1.0}, 
             Table[x[i-1] - 2. x[i] + x[i+1] == 1.0i/n, {i, 2, n-1}], 
             {x[n-1] - 2. x[n] == 1.0}];  
)

You won't want to see all the equations, but here are a few of them.

Short[eqns, 4]

Solve will automatically call the sparse solver.  This example runs many times faster in Version 2.2 than in Version 2.1.

Timing[sol=Solve[eqns, vars];]

Short[sol[[1]],4]

Among the many uses for this functionality is the possibility of 
solving partial differential equations with finite difference methods.

Two point boundary problems: linear problem 
y''=f(x), 0<x<a; y(0)=0=y(a)

Divide interval [0,a] into n subintervals: 0=x(0)<x(1)<...<x(n)=a, and then use finite difference to approximate the derivative y''.  We have 
       y[x(0)]=0; y[x(n)]=0; and  
       y[x(i-1)]-2y[x(i)]+y[x(i+1)]=(a/n)^2 f[x(i)], for 0<i<n

Quit

Set up the equations

n=1000;  
a=5; 
f[x_]=- Sin[x^5]; 
h=N[a/n]; 
vars = Table[y[i], {i, 0, n}]; 
eqns = Join[{y[0] == 0}, 
             Table[y[i-1] - 2. y[i] + y[i+1] == h^2 f[h*i], {i, n-1}], 
             {y[n] == 0}];

sols=Solve[eqns,vars];

ListPlot[Table[{h*i,vars[[i+1]]},{i,0,n}]/.sols[[1]],PlotJoined->True];

Questions: 
y''+ay'+by=f(x), 0<x<a; y(0)=0=y(a);
y''+ay'+by=f(x), with other boundary conditions.

Two point boundary problems: nonlinear problem 
y''=f(x,y), 0<x<a; y(0)=0=y(a)

Divide interval [0,a] into n subintervals: 0=x(0)<x(1)<...<x(n)=a, and then use finite difference to approximate the derivative y''.  We have 
       y[x(0)]=0; y[x(n)]=0; and  
       y[x(i-1)]-2y[x(i)]+y[x(i+1)]=(a/n)^2 f[x(i),y[x(i)]], for 0<i<n

Set up the equations

Quit

n=80;  
a=5; 
f[x_]=- Sin[x^2]; 
h=N[a/n,40]; 
vars = Table[y[i], {i, 0, n}]; 
eqns = Join[{y[0] == 0},Table[y[i-1] - 2 y[i] + y[i+1] == h^2(y[i]^2f[h*i]-1/5), {i, n-1}], 
{y[n] == 0}];

sols=Apply[FindRoot,Join[{eqns},Table[{vars[[i+1]],0.0},{i,0,n}],{WorkingPrecision->20}]];

ListPlot[Table[{h*i,vars[[i+1]]},{i,0,n}]/.sols,PlotJoined->True];

Questions: 
y''=f(x,y,y'), 0<x<a; y(0)=0=y(a);
y''=f(x,y,y'), with other boundary conditions.

Elliptic equation Uxx + Uyy = f(x,y)  
in the rectangle 
       0 < x < a  
       0 < y < b, 
with zero boundary conditions.

Clear[GridU]; 
GridU[n_,m_,a_,b_]:=(dx=N[a/(n+1)]; dy=N[b/(m+1)]; lambda=dx^2/dy^2;  
  Print["Lambda = ",lambda]; 
  Clear[r]; r[i_,j_]=dx^2*f[i*dx,j*dy]; 
  vars = Flatten[Table[u[i,j], {i, n}, {j, m}]]; 
  eqns = Flatten[ 
         Join 
          [ 
           Join[ 
                {-2.0*(1+lambda)*u[1,1] + lambda*u[1,2] + u[2,1] == r[1,1]}, 
                Table[lambda*u[1,j-1] - 2.0(1+lambda)*u[1,j]  
		      + lambda*u[1,j+1] + u[2,j] == r[1,j], {j, 2, m-1}], 
                {lambda*u[1,m-1] - 2.0*(1+lambda)*u[1,m] + u[2,m] == r[1,m]} 
               ], 
           Table[Join[ 
                      {u[i-1,1]-2.0*(1+lambda)*u[i,1] + lambda*u[i,2]  
		       + u[i+1,1] == r[i,1]}, 
                      Table[u[i-1,j]+lambda*u[i,j-1]-2.0*(1+lambda)*u[i,j] 
			    +lambda*u[i,j+1]+u[i+1,j] == r[i,j], {j, 2, m-1}], 
                      {u[i-1,m]+lambda*u[i,m-1] - 2.0*(1+lambda)*u[i,m]  
		       + u[i+1,m] == r[i,m]} 
                     ], 
                 {i,2,n-1}],        
           Join[{u[n-1,1]-2.0*(1+lambda)*u[n,1] + lambda*u[n,2] == r[n,1]}, 
                Table[u[n-1,j]+lambda*u[n,j-1]-2.0*(1+lambda)*u[n,j] 
		      +lambda*u[n,j+1]==r[n,j], {j, 2, m-1}], 
                {u[n-1,m]+lambda*u[n,m-1] - 2.0*(1+lambda)*u[n,m] == r[n,m]} 
               ] 
          ] 
        ] 
    )

nx = 12; ny=24; 
f[x_,y_]=-2x(1-x)-2y(2-y); 
GridU[nx,ny,1.,2.]; 
sols=Solve[eqns,vars]; 
U = (Transpose[Join[ {Table[0.0,{ny+2}]}, 
          Table[Join[{0},Table[u[i,j],{j,ny}],{0}], {i,nx}], 
          {Table[0.0,{ny+2}]} 
        ] /.sols[[1]]]);

Max[U]

ListPlot3D[U,MeshRange->{{0,1},{0,2}}]

tu=Transpose[N[Table[i/13(1-i/13)2j/25(2-2j/25),{i,0,13},{j,0,25}]]];

Dimensions[tu]

Parabolic Equation Ut = c^2 Uxx+f(t,x) 
in	0 < t < T 
	0 < x < a 
with 	U(t, 0) = g1(t) 
	U(t, a) = g2(t) 
	U(0, x) = g(x)

Forward-difference

nt: 	the number of interior points in t-direction 
nx:	the number of interior points in x-direction 
a: 	the length of x-interval 
T:	the time limit

Clear[ParabolicU, g, g1, g2, f]; 
ParabolicU[nt_,nx_,T_,a_,c_]:=( 
    dt=N[T/nt]; dx=N[a/(nx+1)];  
    lambda = c^2*dt/dx^2; Print["Lambda = ",lambda]; 
 
    If[lambda > 0.5, (Print["Lambda > 0.5, the method fails"]; 
			Return[]) 
      ]; 
 
    w[0,0]=g[0]; Do[w[0,i]=g[i*dx],{i,nx}]; w[0,nx+1]=g[a]; 
 
    Do[(w[i,0]=g1[i*dt]; 
	Do[w[i,j]=lambda*w[i-1,j-1]+(1-2*lambda)*w[i-1,j] 
		  +lambda*w[i-1,j+1]+dt*f[i*dt,j*dx], {j,1,nx}]; 
	w[i,nx+1]=g2[i*dt]), {i,1,nt}]; 
 
    U=Table[w[i,j],{i,0,nt},{j,0,nx+1}]; 
    )

g[x_]=Sin[3. x]+x; g1[t_]=0; g2[t_]=1.+t; f[t_, x_]=Sin[x]+t; 
ParabolicU[40,13,1.,N[Pi],1.0]

ListPlot3D[U[[Range[1,40]]],MeshRange->{{0,Pi},{0,1.}}]

Do[ListPlot[U[[i]],PlotRange->{0,3},PlotJoined->True],{i,20}]

Hyperbolic Equation Utt = c^2 Uxx+f(t,x) 
in	0 < t < T 
	0 < x < a 
with 	U(t, 0) = g1(t) 
	U(t, a) = g2(t) 
	U(0, x) = g(x) 
	Ut(0,x) = gt(x)

nt: 	the number of interior points in t-direction 
nx:	the number of interior points in x-direction 
a: 	the length of x-interval 
T:	the time limit

Clear[HyperU, g, gt, g1, g2, f]; 
HyperU[nt_,nx_,T_,a_,c_]:=( 
    dt=N[T/nt]; dx=N[a/(nx+1)]; lambda = N[c^2*dt^2/dx^2];  
    Print["Lambda = ",lambda]; 
 
    If[lambda > 1.,(Print["Lambda is greater than 1, the method fails"]; 
		     Return[]) ]; 
 
    w[0,0]=g[0];  
    Do[w[0,i]=g[i*dx],{i,nx}];  
    w[0,nx+1]=g[a]; 
 
    w[1,0]=w[0,0]+dt*gt[0];  
    Do[(w[1,i]=w[0,i]+dt*gt[i*dx]),{i,nx}];  
    w[1,nx+1]=w[0,nx+1]+dt*gt[a]; 
 
    Do[(w[i,0]=g1[i*dt]; 
	Do[w[i,j]=lambda*w[i-1,j-1]+(2-2*lambda)*w[i-1,j] 
		  +lambda*w[i-1,j+1]-w[i-2,j]+dt^2*f[i*dt,j*dx], {j,nx}]; 
	w[i,nx+1]=g2[i*dt]), {i,2,nt}]; 
 
    U=Table[w[i,j],{i,0,nt},{j,0,nx+1}]; 
    )

g[x_]=Sin[3.14159 x]+1.; gt[x_]=1.; g1[t_]=1.; g2[t_]=1.; f[t_,x_]=(x^4+1)*t^2; 
HyperU[35,15,2.,1.,1.]

Do[ListPlot[U[[i]],PlotRange->{-1,3},PlotJoined->True],{i,35}]