Two examples

One caveat : always try NDSolve first :)

Linear differential equation

LDE=y''[x]+y'[x]-20x y[x]==0;
LDEboundary={y[-1]==5,y[1]==-1};

dsolLDE=DSolve[{LDE,LDEboundary}//Flatten,y,x][[1]]

plot=Plot[{y[x]/.dsolLDE},{x,-1,1},PlotStyle->{Blue,Thick},AxesLabel->{x,y}]

Spectral methods

Define grid and differentiation matrix: [Boyd, Spectral methods, 2000]

initGrid[points_]:=(
nz=points;
z=N[Table[Cos[((r π)/(nz-1))],{r,0,nz-1}],50];
DSpec=Table[If[i==nz&&j==nz,-((2(((nz-1))^2)+1)/6),If[i==1&&j==1,((2(((nz-1))^2)+1)/6),If[i==j,((-z[[i]])/(2(1-((z[[i]])^2)))),((((-1))^(i+j))/(z[[i]]-z[[j]]))((If[i==1||i==nz,2,1])/(If[j==1||j==nz,2,1]))]]],{i,nz},{j,nz}];
D2Spec=DSpec.DSpec;)
initGrid[30]

Check that it works (because mistakes will be there...):

ftest[x_]:=(((x+3/2))^3)Sin[10x](E^(-x));

Show[ListPlot[Thread[{z,ftest[z]}],PlotRange->All,ImageSize->600,PlotStyle->PointSize[Large]],Plot[ftest[x],{x,-1,1},PlotRange->All]]

Check derivatives

fnumd2=D2Spec.ftest[z]//N

Show[ListPlot[Thread[{z,fnumd2}],PlotRange->All,PlotStyle->PointSize[Large]],Plot[ftest''[x],{x,-1,1},PlotRange->All]]

ListPlot[Thread[{z,D2Spec.ftest[z]-ftest''[z]}],PlotRange->All,PlotStyle->PointSize[Large]]

Exponential convergence

Note : now we have to use 50 - digit precise numbers!

initGrid[40]
ListPlot[Thread[{z,D2Spec.ftest[z]-ftest''[z]}],PlotRange->All,PlotStyle->PointSize[Large]]

initGrid[50]
ListPlot[Thread[{z,D2Spec.ftest[z]-ftest''[z]}],PlotRange->All,PlotStyle->PointSize[Large]]

initGrid[30]

Making the matrix

LDE

eqn=D2Spec+DSpec-20 DiagonalMatrix[z];
source=Table[0,{nz}];
eqn//N//MatrixForm;

Naive solution:

LinearSolve[eqn,source]
D2Spec+DSpec-20 DiagonalMatrix[z]//Inverse//N;

Differential equation needs boundary conditions :

LDEboundary

z//N

eqn[[1]]=Table[If[i==1,1,0],{i,nz}];
source[[1]]=-1;
eqn[[-1]]=Table[If[i==nz,1,0],{i,nz}];
source[[-1]]=5;

numsol=LinearSolve[eqn,source];numsol//N
D2Spec+DSpec-20 DiagonalMatrix[z]//Inverse//N;

Check

Show[ListPlot[Thread[{z,numsol}],PlotRange->All,PlotStyle->{Red,PointSize[Medium]}],plot]

ListPlot[Thread[{z,numsol-(y/@z)/.dsolLDE}],PlotRange->All]

Non-linear

DSolve[{y''[x]-y[x]^2-1==0},y,x][[1]]

NDSolve[{y''[x]-y[x]^2-1==0,y[-1]==1,y[1]==2},y,x][[1]]//Timing
plotsol=Plot[y[x]/.%[[2]]/.C[1]->1/.C[2]->1,{x,-1,1}]

Newton’s method

y''[x] - y[x]^2 - 1 = 0

eq:=D2Spec.trial-trial^2-Table[1,{nz}]

trial=1+(z+1)/2;

show:={Show[ListPlot[Thread[{z,trial}]],plotsol],ListPlot[Thread[{z,eq}]],Max[Abs[eq]]//N};
show

jacobian:=(∂_(trial[[i]]))(Subscript[eq, j])

jacobian:=(
matrix=D2Spec-DiagonalMatrix[2 trial];
matrix[[1]]=Table[If[i==1,1,0],{i,nz}];
matrix[[-1]]=Table[If[i==nz,1,0],{i,nz}];
matrix)

trial = trial - LinearSolve[jacobian, {0}~Join~eq[[2 ;; -2]]~Join~{0}]

step:=(
trial=trial-LinearSolve[jacobian,{0}~Join~eq[[2;;-2]]~Join~{0}];
show
)

Table[step,{5}]//Timing;
First[%]
TableForm[%%[[2]]]

Some remarks

Singularities usually have to be extracted analytically

Partial DE not necessarily a problem (30 x 30 matrix -> 900 x 900 matrix, at worst)

Usually coordinate transformation needed (can include ∞)

Periodic problems: use Fourier series