function [yend t y] = shooting(l, q, b) % solves -y'' = l.q(x).y F = @(x,y) [y(2); -l*q(x)*y(1)]; [t y] = ode45(F, [0 b], [0;1], odeset('reltol', 1e-9)); % launch y=0, y'=1 yend = y(end,1); % extract just final value