How to pass pre-generated values to an anonymous function to speed up a fit?
A named function is declared like this:

function out = function_name(in)


If it is written in a file of the same name, and the file is on the MATLAB path, then it is accessible from the base workspace and any function you call. If it is written as a local function (an additional function in a file defining a different function) then it is only accessible to the other functions in that file. If it is in a sub-directory called private, then it is only accessible to the functions defined in the parent directory.

An anonymous function is defined as like this:

function_handle = @(p) p


Here, function_handle is a variable that contains the anonymous function. In other languages this is called a lambda. This lambda can capture variables from its environment:

data = 1:10;
function_handle = @(p) p.*data


function_handle now contains a copy of data. Changing data or clearing it after creating the anonymous function won't affect the copy of data inside it.

The interesting thing for you here is that a named function can also capture data. This is true only for nested functions (which we didn't mention above). The difference between a nested function and a local function is one measly end statement, it's easy to miss it. This is a local function:

function out = main(in)
% code ...
function r = local(p)
% code ...


This is identical to:

function out = main(in)
% code ...
end
function r = local(p)
% code ...
end


The end statements ending a function are optional.

This is a nested function:

function out = main(in)
% code ...
   function r = nested(p)
   % code ...
   end
end


Here the end statements are mandatory. nested is inside main. nested also has access to variables declared in main. For example, nested can use in.

main can also create a function handle of nested and return it. This function handle then will point at the nested function, and capture its environment. Much like an anonymous function:

function out = main(in)
   data = integral(in,0,1);
   out = @r;
   function r = nested(p)
      r = p .* data;
   end
end


Now, calling main returns a handle to the function nested, which has captured data computed in main. You can use this to create a function handle that you can pass to fit. nested would be your zeusEquation, and main would simply package it with the precomputed integral.

Code review

There are a few things to note about your code:


You compute sin(psi0) and 1-cos(psi0) repeatedly inside your loop. psi0 is a constant within this loop, and consequently you should compute these values once outside of the loop.
You compute integral(sinIntegral,0,u(a6,a5,x(i))) twice for each x(i). This being an expensive function, you should compute it only once and re-use the result.
for i = 1:max(size(x)) is long-hand for for i = 1:length(x). This only works correctly in your loop if x is a vector. In this case, you can use numel(x) instead of length(x). numel is faster.


These three points above lead to:

sinpsi0 = sin(psi0);
cospsi0_1 = 1-cos(psi0);
for i = 1:numel(x)
   if u(a6,a5,x(i))~=0
      int = integral(sinIntegral,0,u(a6,a5,x(i)));
      y(i) = a3 + a4.*(x(i)-a5)  ...
                + (1-int) .* (a1.*sign(x(i)-a5).*sinpsi0 ...
                              - (1-a2).*int.*cospsi0);
   else
      y(i) = a3 + a4.*(x(i)-a5);
   end
end



integral must be called with scalar limits, but other components can be computed vectorised. It would be possible to compute only the integral within the loop, and the rest without a loop.

This leads to:

int = zeros(size(x));
for i = find(u(a6,a5,x)~=0)
   int(i) = integral(sinIntegral,0,u(a6,a5,x(i)));
end
y = a3 + a4.*(x-a5)  ...
       + (1-int) .* (a1.*sign(x-a5).*sin(psi0) ...
                     - (1-a2).*int.*(1-cos(psi0)));


This should be faster, but I haven't tested it. (I assume a1 and the like are all scalars?)
You define your functionzeusEquation using varargin, then proceed to unpack the varargin cell array with the assumption that all your input arguments are there. This is correct in your situation of course, but you might as well do:

function y = zeusEquation(psi0,a1,a2,a3,a4,a5,a6,x)


In this way, if zeusEquation is called with fewer (or more!) arguments, you'll get a sensical error message, rather than an "index out of bounds" error.