Sorry Jack. What do you mean? On my Computer there are no problems.

Regards Hans

Attachments:

Seems your function does not depend too much on r.....

Clear[ff]
ff[r_] := sol[1, r, 3][[1, 1]] /. Rule[x_, y_] :> y[xx]

Manipulate[
 Plot[ff[r], {xx, 0, 1}, PlotRange -> {0, 2}],
 {{r, 0}, -5, 5}]

In[119]:= (f /. sol[1, -5, 3][[1, 1]])[.546]
(f /. sol[1, -1, 3][[1, 1]])[.546]
(f /. sol[1, 0, 3][[1, 1]])[.546]
(f /. sol[1, 2.3, 3][[1, 1]])[.546]
(f /. sol[1, 4, 3][[1, 1]])[.546]

Out[119]= 1.04275

Out[120]= 1.04316

Out[121]= 1.04326

Out[122]= 1.0435

Out[123]= 1.04367

some mistypings.....

eqns[a_?NumericQ, r_?NumericQ, 
  gamma_?NumericQ] := {Derivative[4][f][y] - 
    2 Derivative[2][f][y] a^2 + a^2 (-r y + a^2) f[y] == 0, f[0] == 0,
   Derivative[3][f][0] - Derivative[1][f][0] a^2 == 0, 
  Derivative[1][f][0] == 1, Derivative[2][f][0] == gamma}


sol[a_?NumericQ, r_?NumericQ, \[Gamma]_?NumericQ] := 
 NDSolve[eqns[a, r, \[Gamma]], f, {y, 0, 1}, Method -> "BDF", 
  MaxSteps -> Infinity, InterpolationOrder -> All]

eval = sol[1, 2, 3][[1, 1]]

ff = f /. eval

Plot[ff[y], {y, 0, 1}]