# Beginners problem, Do Loop, Eigenfunction iteration

I am trying to find the first Eigenfunction of the Laplacian (in 1D), i.e. a solution of $$u''(x)=k u(x)\\ u(0)=u(1)=0$$ with minimal $k>0$ (in this trivial example, I actually know the analytic solution but this is not the point; neither is using any built-in Eigensolver. I simply want to get acquainted with Mathematica).

To do so, I want to use the Power iteration method:

Repeatedly applying the inverse Laplacian to an arbitrary initial function will produce a sequence converging to the solution.

What I currently have is

ClearAll[h];
ClearAll[f];
f[x_] = 0.5 - Abs[x - 0.5];
steps = 1;
Do[
s = NDSolve[{h''[x] == f[x], h[0] == 0, h[1] == 0}, h, {x, 0, 1}];
f[x_] = h[x]/h[0.5] /. s;
, steps]

However, when I take steps=2 or more, I get

Dot::rect: Nonrectangular tensor encountered.

Note that this is my very first use of Mathematica, and you might have to explain obvious things to me.

Side question: How do I properly output intermediate results in the loop? If I use

f[x_] = h[x]/h[0.5] /. s; ?f

I do get info on f but if I put

f[x_] = h[x]/h[0.5] /. s; ?f ; ?s

I get

Information::nomatch: "No symbol matching ?s found."

(Adding newlines or using commas doesn't help either. For example,

f[x_] = h[x]/h[0.5] /. s;
?f ,
?s ,

gives no output at all, whereas

f[x_] = h[x]/h[0.5] /. s;
?f ;
?s ;

gives

Information::ssym: "\!$$f; Information[\"s\", LongForm -> False]$$ is not a symbol or a valid string pattern"

)

• So at each step the definition of f and s will be interpolating functions. What information do you want to have about them? How do you want to display it? Commented May 11, 2016 at 10:43

To solve the first issue, write

s = NDSolve[{h''[x] == f[x], h[0] == 0, h[1] == 0}, h, {x, 0, 1}][[1]];

instead (notice the final brackets).

• …or, use First[]. Commented May 11, 2016 at 10:39

This allows you to visualize the (small) changes to the function over the course of the iterations,

ClearAll[h];
ClearAll[f];
f[x_] = 0.5 - Abs[x - 0.5];
steps = 4;
list = Join[
{{#, f[#]} & /@ Subdivide[1, 100]},
Table[
s = NDSolveValue[{h''[x] == f[x], h[0] == 0, h[1] == 0},
h, {x, 0, 1}];
f[x_] = s[x]/s[0.5];
{#, f[#]} & /@ Subdivide[1, 100]
, steps]];
ListLinePlot[list, PlotLegends -> Range[0, 4]]

I've also used NDSolveValue to get around your first problem.