I found the solution:
syms x1 x2 z1 z2 a u
% Defining x_dot = f(x)
f(x1,x2) = [-2*x1+a*x2+sin(x1);
-x2*cos(x1)+u*cos(2*x1)];
% Defining z = g(x)
g(x1,x2) = [x1;
a*x2+sin(x1)];
solver = formula(g);
% Defining J as the jacobian of g: J = dg(x)/dt = dz/dt
J(x1,x2) = jacobian(g);
% Defining x as a function of z: x = g^{-1}(z)
[x1, x2] = solve(z1 == solver(1), ...
z2 == solver(2));
% Calculate z_dot = dg(x)/dx * dx/dt | x = g^{-1}(z)
% Note: dg(x)/dt = dz/dt
z_dot = J(x1,x2)*f(x1,x2);
simplify(z_dot)
which outputs:
ans =
z2 - 2*z1
cos(z1)*sin(z1) - 2*z1*cos(z1) + a*u*(2*cos(z1)^2 - 1)
And since:
>> simplify(2*cos(z1)^2 - 1)
ans =
cos(2*z1)
Then the final solution is:
ans =
z2 - 2*z1
cos(z1)*sin(z1) - 2*z1*cos(z1) + a*u*cos(2*z1))
Which is the correct solution.
