how to code for system of simultaneous nonlinear differential equation?

xt^2 dy/dt +yt^2 dx/dt +2xyt=0;

x dx/dt +y^k dy/dt +GM/t^2 =0;

all the derivatives are partial, since x=x(t,r),y=y(t,r);

r is not used here but x and y are dependent on it.

k lies between 1 and 5/3;

G is gravitational const.

M is mass of any massive object,say sun.

[Walter Roberson]

x=x(t,r)

r is not used here but x and y are dependent on it.

It is not clear to me whether you are working with a Partial Differential Equation system or not? Would there be plans to vary r in the future, or will you intentially only be computing along lines of constant r at any given time ?

[Walter Roberson]

I do not seem to have usable initial conditions.

You cannot start at time t = 0 because there is a division by time implied

syms GM k positive
assume(k>=1 & k<=5/3)
syms x(t) y(t)
dx = diff(x); dy = diff(y);

eqn1 = x * t^2 * dy + y * t^2 * dx + 2 * x * y * t == 0;
eqn2 = x * dx + y^k * dy + GM/t^2 == 0;
eqns = [eqn1, eqn2]

eqns(t) = 

these_eqns = subs(eqns, {GM, k}, {1e20, 5/4})

these_eqns(t) = 

[eqs,vars] = reduceDifferentialOrder(these_eqns,[x(t), y(t)])

eqs = 

vars = 

initConditions = [1e5 2e7];
[M, F] = massMatrixForm(eqs, vars)

M = 

F = 

f = M\F

f = 

odefun = odeFunction(f,vars)

odefun = function_handle with value:
    @(t,in2)[(1.0./t.^2.*in2(1,:).*(t.*in2(2,:).^(9.0./4.0)-5.0e+19).*2.0)./(in2(1,:).^2-in2(2,:).^(9.0./4.0));(1.0./t.^2.*in2(2,:).*(t.*in2(1,:).^2-5.0e+19).*-2.0)./(in2(1,:).^2-in2(2,:).^(9.0./4.0))]

[t, out] = ode15s(odefun, [1 100], initConditions);

Warning: Failure at t=1.000118e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

plot(t, out)
legend({'x', 'y'})