function A6_makeall
  clear all;

  %define the interval, we want to see
  X = -1:0.1:1;
  Y = -1:0.1:1;

  %calc the Lyapunov Function
  for myi = 1:size(X, 2)
    for myj = 1:size(Y, 2)
      Lyapunov_Function(myi, myj) = (X(myi)^2 + Y(myj)^2)/2;
    end
  end
  
  %solve the system
  [T, X_solution] =  ode45(@diffx, ...    % function for the system
                           [0, 1000], ... % time
                           [1, 0]);       % initial conditions

  %calc the value of the Lyapunov Function along the solution
  for myi = 1:size(X_solution(:, 1))
    Lyapunov_Function_along_solution(myi) ...
        = (X_solution(myi, 1)^2 + X_solution(myi, 2)^2)/2;
  end

  meshc(X, Y, Lyapunov_Function);           
  hold all;                       % plot everything in one figure
  plot(X_solution(:, 1), X_solution(:, 2)); 
  plot3(X_solution(:, 1), X_solution(:, 2), ...
        Lyapunov_Function_along_solution);  

end


function xdot = diffx(t, x)

  xdot_1 = -x(1)^3 + 2*x(1)*x(2);
  xdot_2 = -x(1)^2 - x(2);

  xdot = [xdot_1; xdot_2];
end