%% Równanie diofantyczne
am0=Am(1)
am1=Am(2)
am2=Am(3)
am3=Am(4) % elementy Am (mianownik z nowymi biegunami)
a0=Den(1)
a1=Den(2)
a2=Den(3)
a3=Den(4) % elementy A (mianownik oryginalny)
b0=Num(1)
b1=Num(2)
b2=Num(3)
b3=Num(4) % elementy B (licznik oryginalny)


eqns = [
 am1 == f1*a0 + a1
 am2 == f2*a0 + f1*a1 + a2
 am3 == f3*a0 + f2*a1 + f1*a2 + a3
 0 == f4*a0 + f3*a1 + f2*a2 + f1*a3
 0 == f5*a0 + f4*a1 + f3*a2 + f2*a3
 0 == f6*a0 + f5*a1 + f4*a2 + f3*a3
 0 == f7*a0 + f6*a1 + f5*a2 + f4*a3
 0 == f8*a0 + f7*a1 + f6*a2 + f5*a3
 0 == f9*a0 + f8*a1 + f7*a2 + f6*a3
 0 == f10*a0 + f9*a1 + f8*a2 + f7*a3
 0 == f11*a0 + f10*a1 + f9*a2 + f8*a3
 0 == f12*a0 + f11*a1 + f10*a2 + f9*a3
 0 == f13*a0 + f12*a1 + f11*a2 + f10*a3
 0 == f14*a0 + f13*a1 + f12*a2 + f11*a3 + b0*g0
 0 == f15*a0 + f14*a1 + f13*a2 + f12*a3 + b1*g0 + b0*g1
 0 == f16*a0 + f15*a1 + f14*a2 + f13*a3 + b2*g0 + b1*g1 + b0*g2
 0 == f17*a0 + f16*a1 + f15*a2 + f14*a3 + b3*g0 + b2*g1 + b1*g2
 0 == f17*a1 + f16*a2 + f16*a3 + b3*g1 + b2*g2
 0 == f17*a2 + f16*a1
 0 == f17*a3
];


R=solve(eqns)
F(1) = 1;
for i = 1:dF
 F(i+1) = eval(['double(R.f',num2str(i),')']);
end
for i = 0:dG
 G(i+1) = eval(['double(R.g',num2str(i),')']);
end