Since you have not given the set of differential equations, this answer cannot respond exactly. However, as indicated in the comments, the system of equations is implicit. These can be written as $F(t,y,y')=0$, where $$y(t)=[y_1(t),y_2(t),\ldots,y_n(t)]^T \text{ and } y'(t)=[y'_1(t),y'_2(t),\ldots,y'_n(t)]^T.$$This system can be solved numerically in MATLAB using ode$15$i.
As an example, suppose we want to solve the implicit system of ODEs
$$
F(t,y,y')=\begin{bmatrix}\cos(y'_1)+y'_2-y_1\\t\sin(y'_2)+y'_1-y_2^2\end{bmatrix}
$$
with $y(0)=[0,0]^T$ for $0\leq t\leq T$.
ode$15$i requires consistent initial conditions for both $y$ and $y'$. Given that we do not have initial conditions, they can be approximated using the decic command.
This can be done in MATLAB as follows.
T = 5;
y0 = [0;0];
fun = @(t,y,yp) [cos(yp(1))+yp(2)-y(1);t*sin(yp(2))+yp(1)-y(2)^2];
[y0,yp0] = decic(fun,0,y0,[1;1],[0;0],[0;0]);%get consistent initial conditions
[tsol,ysol] = ode15i(fun,[0,T],y0,yp0); % solve the system
plot(tsol,ysol(:,1),'b',tsol,ysol(:,2),'r') % plot the solution
Edit: Solution to include equations
System
The system of DEs to be solved is as follows:
$$
\begin{align}
\theta '' &= \frac{\theta' \cos \theta}{R} + \frac{R'\sin \theta}{R} + \frac{\nu\sin\theta}{\kappa R} + \frac{\eta\cos\theta}{\kappa R} \tag{1}\\
\nu' &= \frac{\kappa}{2} \left( \theta' + \frac{\sin\theta}{R}-c_0 \right)\left( \theta'-\frac{\sin\theta}{R}-c_0\right) + \lambda\tag{2}\\
\eta' & =0\tag{3}\\
R' &= \cos \theta \tag{4}\\
z' &=- \sin \theta \tag{5}\end{align}
$$
No initial conditions are given, though it is stated that all zero initial conditions may be appropriate.
Solution
Noting $(3)$, the solution of $\eta'=0$ must be $\eta(t)=c$ for some constant $c$. Since $\eta(0)=0$, we must have $c=0$, giving $\eta(t)=0$. As such, $(3)$ is not needed, since $\eta(t)$ can be obtained analytically. I will assume, however, throughout the remainder of the solution below, that $\eta(t)$ need not be a zero constant, but simply any constant, as given by the initial conditions.
Furthermore, $(1)$ can be reduced to first order by introducing the auxiliary variable $A=\theta'$. Then $(1)$ becomes
$$
A'=\frac{A \cos \theta}{R} + \frac{R'\sin \theta}{R} + \frac{\nu\sin\theta}{\kappa R} + \frac{\eta\cos\theta}{\kappa R} \tag{6}
$$
with the additional equation
$$
\theta'=A\tag{7}.
$$
The system then required to be solved is $(2),\,(4),\,(5),\,(6)\text{ and }(7)$.
In MATLAB, an m-file for the ode function should be created similar to as follows
function F = odefun(t,y,yp,parms)
% evaluate the ode functions with the given parameters
% y and yp are ordered as [theta,nu,R,z,A]
% parms is a structure, containing eta, kappa, lambda, c0
% extract components from parms structure
eta = parms.eta;
kappa = parms.kappa;
lambda = parms.lambda;
c0 = parms.c0;
% extract components from y and yp for more clear use
theta = y(1); thetap = yp(1);
nu = y(2); nup = yp(2);
R = y(3); Rp = yp(3);
z = y(4); zp = yp(4);
A = y(5); Ap = yp(5);
% allocate memory for return variable
F = zeros(5,1);
% precompute some quantities
sint = sin(theta);
cost = cos(theta);
% evaluate equation (2)
F(1) = nup - (kappa/2*(thetap+sint/R-c0)*(thetap-sint/R-c0)+lambda); %****
% evaluate equation (4)
F(2) = Rp - cost;
% evaluate equation (5)
F(3) = zp + sint;
% evaluate equation (6)
F(4) = Ap - (A*cost + Rp*sint + (nu*sint+eta*cost)/kappa)/R; %****
% evaluate equation (7)
F(5) = thetap - A;
Note: Both equations marked by %**** $\left(\text{i.e. $(2)$ and $(6)$}\right)$ will have a divide by zero if $R(0)=0$.
From my previous smaller example, you should be able to make this function work with the ode15i(...) and decic(...) commands.
To avoid the need of algebraic manipulation with $s,t$, you can use the folowing method. We have:
$$\frac{-dp}{2pz} = \frac{-dq}{2qz} \implies \frac{dp}{p} = \frac{dq}{q} \implies p = Cq$$
Now plugging into the equation we have that $C^2q^2 + q^2 = z^2 \implies q = \frac{z}{\sqrt{C^2+1}}$ and $p = \frac{Cz}{\sqrt{C^2+1}}$. Then plugging into $dz = pdx + qdy$ and integrating we have:
$$dz = \frac{Cz}{\sqrt{C^2+1}}dx + \frac{z}{\sqrt{C^2+1}}dy$$
$$\frac{dz}{z} = \frac{C}{\sqrt{C^2+1}}dx + \frac{1}{\sqrt{C^2+1}}dy$$
$$\ln|z| = \frac{Cx + y}{\sqrt{C^2+1}} + D$$
$$z = Ae^{\frac{Cx + y}{\sqrt{C^2+1}}}$$
This is the complete integral of the PDE.
Now for the particualr solution we use:
$$\frac{d}{dt}(1) = u_x\frac{dx}{dt} + u_y\frac{dy}{dt}$$
$$0 = C(-\sin t) + \cos t \implies C(t) = \cot t$$
Plugging this with the initial condition in the complete integral after some manipulation we have that
$$1 = A(t)e \implies A(t) = e^{-1}$$
Finally we get $u(x,y) = e^{\cos(t)x + \sin(t)y -1}$, where $t \in [0,2\pi]$. This is a family of solution that passes through the curve. To find a single solution that passes through the curve we need to find the envelope of the family of solutions given by $\phi(x,y,u,A(t),C(t)) = 0$. This happens when $\frac{d\phi}{dt} = 0$ and we get:
$$0 = \frac{d\phi}{dt} = \phi_A \frac{dA}{dt} + \phi_C \frac{dC}{dt} = \frac{x - Cy}{(1+C^2)^{\frac 32}} \implies x = yC$$
Plugging into the family of soltions we get the wanted solution:
$$\boxed{u(x,y) = e^{\sqrt{x^2+y^2} - 1}}$$
REMARK: From the solution we can conclude that if we write it in polar coordinates it will be independent of $\theta$. This prompts you to believe that if we introduce change of variables at the beginning then we would have an easier job. As an exercises maybe can check this.
Best Answer
As I have pointed out in my comment, this system can be reduced to a more treatable one. Once you have cast it into the form
$$\left\{ \begin{array}{l} \xi^\prime=\xi^2+b\xi+c\eta+d\theta\\ \eta^\prime=\eta^2+f\xi+g\eta+h\theta\\ \theta^\prime=\theta^2+j\xi+k\eta+l\theta \end{array} \right.$$
one introduces the three new variables
$$\left\{ \begin{array}{l} \xi=-\frac{1}{X}\frac{dX}{dt} \\ \eta=-\frac{1}{Y}\frac{dY}{dt} \\ \theta=-\frac{1}{Z}\frac{dZ}{dt} \end{array} \right.$$
and you will remove the quadratic terms, obtaining
$$\left\{ \begin{array}{l} \frac{1}{X}\frac{d^2X}{dt^2}=b\frac{1}{X}\frac{dX}{dt}+c\frac{1}{Y}\frac{dY}{dt}+d\frac{1}{Z}\frac{dZ}{dt} \\ \frac{1}{Y}\frac{d^2Y}{dt^2}=f\frac{1}{X}\frac{dX}{dt}+g\frac{1}{Y}\frac{dY}{dt}+h\frac{1}{Z}\frac{dZ}{dt} \\ \frac{1}{Z}\frac{d^2Z}{dt^2}=j\frac{1}{X}\frac{dX}{dt}+k\frac{1}{Y}\frac{dY}{dt}+l\frac{1}{Z}\frac{dZ}{dt} \end{array} \right.$$
This system can be easily rewritten as
$$\left\{ \begin{array}{l} \frac{d^2X}{dt^2}=b\frac{dX}{dt}+X\frac{d}{dt}\ln(Y^cZ^d) \\ \frac{d^2Y}{dt^2}=g\frac{dY}{dt}+Y\frac{d}{dt}\ln(X^fZ^h) \\ \frac{d^2Z}{dt^2}=l\frac{dZ}{dt}+Z\frac{d}{dt}\ln(X^jY^k) \end{array} \right.$$
that can be solved with standard exponentials provided a relation between the coefficients exists so that the derivatives of the logarithms are zero.
For the more general case I cannot be of help but one can always reduce to a similar set of equations.