Tak jsem si trochu hrál:
clc
clear all
% ---------------
krok=100;

%reseni
reseni = solve('y=x+1','y=x^2-2*x-1');
prusecik_x=double(reseni.x) %musi se konvertovat na double, jinak by to byl "sym object" a s tim se nada moc pracovat
prusecik_y=double(reseni.y)

x=round(prusecik_x(2)):1/krok:round(prusecik_x(1)); %generuje x-ove hodnoty, nejprve si to zaokrouhlim, aby to byly pekny cisla
y=x.^2-2*x-1; %vypocet bodu paraboly-bacha na (.^) ta tecka urcuje ze jde o Array power
y2=x+1; %vypocet bodu primky

% vykresleni:
plot(x,y,... paraboly
x,y2,... primky
prusecik_x(1), prusecik_y(1),'or',... prvniho pruseciku - "o" znaci, ze to bude kolecko a "r" jako red
prusecik_x(2), prusecik_y(2),'or') %druheho pruseciku

title('zobrazeni pruseciku primky a paraboly danymi rovnicemi y=x+1, y=x^2-2x-1') %nazev grafu
xlabel('x-ova osa') %popisek x-ove osy
ylabel('y-ova osa') %popisek y-ove osy
grid %mrizka

%Pro kontrolu analyticky vypocet:
%Kdyz dosadis za "y" do rce y=x^2-2*x-1 tu rci primky y=x+1 vyleze ti
%rce:x^2-3x-2
disp('------------------------')
disp('analyticke reseni:')
disp('------------------------')
x_koreny=roots([1 -3 -2])
y_koreny=x_koreny.^2-2*x_koreny-1
Stačí to okopčit a vložit do čistého m-file, snad je to to, co jsi chtěl. Celkem jednoduše bys to upravil na řešení libovolných parabol a přímek.