Thank you Doug and Robert! I will try to follow the paths shown. It really helps.

Hi there, I tried all these methods mentioned above before I pasted. The methods are only valid when f(x1,x2,x3,x4) is analytical and independent. The problem is that the f in my code does not exist by itself. It comes from the numerical solution which actually treat x2,x3,x4 as parameters. This is the source of the error message. x0:= 1/2: s:=x->x^2: sol:=(x2,x3,x4)->dsolve({y''(x1)-(x3^2+s(x1)*x4^2)*y(x1)=0,y(x2)=exp(-x2),y'(x2)=1/2},numeric,output=operator,range=0..x0): g:=(x1,x2,x3,x4)->eval(y(x1),sol(x2,x3,x4)): dx2:=1/1000: Data:=(x1,x2,x3,x4)->[seq([x2+i*dx2,g(x1,x2+i*dx2,x3,x4)],i=-3..3)]: f := (x1, x2, x3, x4) -> PolynomialInterpolation(Data(x1, x2, x3, x4), x2) : dx2Data := (x1, x2, x3, x4) -> eval(diff(f(x1, x, x3, x4), x), x = x2) dx2Data((1/3)*x0, (1/3)*x0, 10, 1); Error, (in dsolve/numeric/process_input) input system must be an ODE system, got independent variables {x1, x-3/1000} Mathematica can do this smoothly. But I can not get the right answer in Maple. Can anyone help? I am desperate about this. Cheers.