Is there a way to solve a system of nonlinear equations without the optimization toolbox (fsolve) or symbolic math toolbox (sym)?

I do not have access to Optimization or Symbolic Math Toolboxes but I need to solve a system of 2 nonlinear equations with 2 unknowns. All of the examples I have come across so far have included some combination of "fsolve", "solve", "syms", "root2d", etc and I do not have access to any of those. Is there any way to solve this system without using those 2 toolboxes? My equations themselves are below and the 2 unknown variables should only have 1 solution. I am also running Matlab R2014b if that helps

% Simplifying equation (all known constants)
A = 4*kappa
B = 2*alpha/(rho*c*A)
D = -(r^2)

% Equations
T1 = B*x*y(-expint(-D/x)+expint(D/(x+A*t1)))
T2 = B*x*y(-expint(-D/x)+expint(D/(x+A*t2)))

% A, B, D, T1, T2, t1, t2 are all known
% x and y are unknown 

Answer 1

[Walter Roberson]

Systems of equations can often be turned into a minimization.

A = 4*kappa
B = 2*alpha/(rho*c*A)
D = -(r^2)
% Equations
eqn1 = @(x,y)T1 - B.*x.*y.*(-expint(-D./x)+expint(D./(x+A*t1)))
eqn2 = @(x,y)T2 - B.*x.*y(-expint(-D./x)+expint(D./(x+A*t2)))
residue = @(xy) eqn1(xy(1),xy(2)).^2 + eqn2(xy(1),xy(2)).^2;

xy0 = rand(1,2);  %some initial guess
[bestxy, fval] = fminsearch(residue, xy0)
bestx = bestxy(1); besty = bestxy(2);

At this point, fval would ideally be 0. If it is substantially different than 0 then the solution from fminsearch is not good and you might want to try a different xy0 .