Hi folks, I am solving a boundary value problem numerically using dsolve. I would like to pass the resulting solution (which I gather is really a procedure) back into another differential equation as a coefficient.

M:=x(xi)^2+1:

ode1:=diff(M*diff(x(xi),xi),xi)-0*x(xi)=0:

ics:=x(0)=0,x(1)=1:

sol1:=dsolve({ode1,ics},numeric,output=listprocedure):

y_sol1:=sol1[2]:

y_sol1:=rhs(y_sol1):

% I could then plot the solution using