function [ub, uu, x]=u_ivp(lambda, a, b, ua, reltol, abstol, Npts) IVs = [ua; 1]; opt = odeset('RelTol',reltol,'AbsTol',abstol); [x u] = ode45(@rhs, linspace(a,b,Npts), IVs, opt, lambda); % u(b) and u(x) ub = u(end,1); uu = u(:,1);