Trouble solving system of equations -- trying vpasolve

I would like to solve for parameters a, b, c, and d in the following function:

y(x, a, b, c, d) = a*exp(b*x + c) + d

I get a system of four equations by constraining the values of the function and its derivative at two points.

y1 = y(x1, a, b, c, d)

dydx1 = dydx(x1, a, b, c, d)

y2 = y(x2, a, b, c, d)

dydx2 = dydx(x2, a, b, c, d)

I am trying to use vpasolve in the code given below. To test, I choose values for the parameters and for x1 and x2 and calculate y1, y2, dydx2, and dydx2. I use x1, x2, and the calculated values to set up the system of equations. If I run vpasolve with no initial guess, I get an answer I don't understand. If I use the known solution as the initial guess, I get no answer.

Do I have a problem with my script? Should I do something different? I get the result I expect in a simpler example with only two parameters. In the simpler case, I set c and d to zero and only constrain by the function values and not the derivatives.

To clarify, in the code given below, I hoped vpasolve would find the known solution to the system of equations as defined in the script, aa=aaa, bb=bbb, cc=ccc, dd=ddd. This doesn't happen even if I provide the known solution as the initial guess. I'm asking for help with figuring out what is going on and for advice about what I could try to get better results.

I'm using MATLAB R2019a.

Thank you!

clear variables;

%% define symbolic variables and equations
syms y(x, a, b, c, d) aa bb cc dd 

% system of equations will constrain using function and derivative
y(x, a, b, c, d) = a*exp(b*x + c) + d;
dydx = diff(y, x);

%% create test case with known solution
aaa = 1;
bbb = -1;
ccc = 1;
ddd = 0;

% get values of function and derivative at two points
x1 = 0;
y1 =    eval(   y(x1, aaa, bbb, ccc, ddd));
dydx1 = eval(dydx(x1, aaa, bbb, ccc, ddd));

x2 = 1;
y2 =    eval(   y(x2, aaa, bbb, ccc, ddd));
dydx2 = eval(dydx(x2, aaa, bbb, ccc, ddd));

%% set up system of equations and solve
eq1 =    y(x1, aa, bb, cc, dd) == y1;
eq2 = dydx(x1, aa, bb, cc, dd) == dydx1;
eq3 =    y(x2, aa, bb, cc, dd) == y2;
eq4 = dydx(x2, aa, bb, cc, dd) == dydx2;

eqns = [eq1, eq2, eq3, eq4].';
vars = [aa, bb, cc, dd].';

% try setting initial guess to be known solution
initial_guess = [aaa, bbb, ccc, ddd].';
%initial_guess = [];

sol = vpasolve(eqns, vars, initial_guess);  

%% test solution
isAllGood = all(isAlways(subs(eqns, sol)));

%% plot solution
xx = x1:0.1:x2;
if ~isempty(sol.aa)
    yy = eval(y(xx, sol.aa, sol.bb, sol.cc, sol.dd));
else
    yy = NaN(size(xx));
end
yyTrue = eval(y(xx, aaa, bbb, ccc, ddd));

figure(1);
clf;
hold on;
plot(xx, yy);
plot(xx, yyTrue, '--');
plot(x1, y1, 'xr');
plot(x2, y2, 'xr');
grid on;
xlim(sort([x1 x2]) + [-1 1]);
ylim(sort([y1 y2]) + [-1 1]);
legend({'From solver' 'Truth'});
xlabel('x');
ylabel('f(x)');

Answer 1

[darova]

My attempt using fsolve. Result depend on initial guess

clc,clear

syms a b c d x
p = [a b c d];
% system of equations will constrain using function and derivative
f = p(1)*exp(p(2)*x + p(3)) + p(4);
df = diff(f,x);
F = matlabFunction(f);
F = @(p,x) F(p(1),p(2),p(3),p(4),x);
dF = matlabFunction(df);
dF = @(p,x) dF(p(1),p(2),p(3),x);

%% create test case with known solution
a1 = 1;
b1 = -1;
c1 = 1;
d1 = 0;
p1 = [a1 b1 c1 d1];
x = 0:0.1:1;
y = F(p1,x);
plot(x,y)

eqs = @(x)[F(x,0)-F(p1,0)
            F(x,1)-F(p1,1)
            dF(x,0)-dF(p1,0)
            dF(x,1)-dF(p1,1)];
opt = optimset('MaxFunEvals',3000,'MaxIter',1000);
p2 = fsolve(eqs,ones(1,4)*3,opt);

hold on
plot(x,F(p2,x),'r')
hold off