There is nothing wrong with the approach that you want to use, however there are some details that will help improve the efficiency of the method.
First, I the really smart Matlab people are always reminding me that the variable i is a defined Matlab constant, and is used in complex math notation, therefore, using i as a loop index can be considered a bad habit. Your code will work using i as the counting variable, or loop index, but just realize that you are ghosting a Matlab constant.
Next, it is my opinion that a counting loop should always be structured such that when the loop ends, the counting variable contains the number of loop itterations (and the number of elements in the accumulated vectors). To accomplish this, the counter should be incremented at the start of the loop, and therefore it needs to be initialized to zero outside the loop; i.e.:
i=0
while (..)
i=i+1; % this is the first statement inside the loop
...
...
end
Now when the loop is completed, i contains the number of elements in x1, y1, and t1.
Next, your algorithm will be very sensitive to the step size because you are using the last loop itteration as the final value for the distance. This means that there will be an error in the final value that is a function of the step size, and the only way to minimize this error is to make the step very small. A better aproach is to use linear interpolation between the last and the next-to-last step to fine the precise point where the vertical trajectory crossed the zero boundary (i.e. the ground plane). Using interpolation will allow you to get a more precise solution while using a much larger step size.
To do this, you need to change the termination condition on the while lop so that it will always terminate on the first step after the Y position passes zero (goes negative). You can do this using
while y1(i) >= 0
end
(edit: on the first pass, y1(i) will have a subscript of zero, and this will give an error. Need to change this to avoid the error)
h = 0
while h>=0
i=i+1;
...
...
h=y1(i);
end
This will allow you to initialize y1 at zero and the loop will not terminate until y1 goes negative.
After the loop terminates, you want to replace the last value in x1,y1,and t1 with the interpolated value.
First, calculate the interpolation fraction (remember i contains the number of points in the trajectory, and the last point is below the ground):
k = -y1(i-1) / (y1(i) - y1(i-1));
Now calculate the interpolated values:
xinterp = x(i-1) + k * (x(i) - x(i-1));
yinterp = 0;
tinterp = t(i-1) + k * (t(i) - t(i-1));
Now replace the end values with the iterpolated ones:
x1(i) = xinterp;
y1(i) = yinterp;
t1(i) = tinterp;
Now the last point in the trajectory is exactly on the ground.
the length of the jump = x1(i)
the time of the jump = t1(i)
and the height of the jump = max(y1)
After making this change, you should be able to get a pretty good answer using a much biger step size, and this will speed up your code.
Finally, in order to locate the initial velocity for a given distance, use an interpolation method to calculate the next guess and you will converge on the answer much faster than simply taking small steps in the initial velocity.
