Alireza,
As you've noted, R gets overwritten on each pass of the loop. Unfortunately I don't have the symbolic math toolbox so I can't test your code exactly, but the indexing should work since x is a scalar. In the code you posted, you're assigning the initial guess to R instead of the solved value. Perhaps this is the issue?
Try this. I've preallocated R prior to the loop (a bit more efficient). It should give you a column vector of symbolic variables matching the (row) vector of xinitial.
syms x
y = 1 + cos(x) * cosh(x)
xinitial = linspace(0,30,61);
n = length(xinitial);
R = sym('0',[n,2])
for i = 1 : n
R(i,1) = vpasolve(y == 0, x, xinitial(i))
end
