Solve boundary value problem y''+y = f, 0<x<1, y(0)=y(1)=0

restart:

with(plots):

Find sine series coefficients of f(x)

assume(n,integer):

f:=x->x^3-x;

N:=10:

for n from 1 to N do

b[n]:=2*int(f(x)*sin(n*Pi*x),x=0..1);

od:

Calculate Fourier sine series of solution

y:=x->add(b[n]/(1-n^2*Pi^2)*sin(n*Pi*x),n=0..N):

plot(y,0..1);

Find exact solution

exact:=rhs(dsolve({diff(z(x),x$2)+z(x)-f(x),z(0)=0,z(1)=0},z(x)));

Plot error

plot(y(x)-exact,x=0..1);