Fract := proc(P::posint, Q::posint)  
         local p,q:
         for p from 1 to P-1 do
            for q from 1 to Q-1 do
              if (P-p)*q-P*(Q-q)=1 the return (p/q,(P-p)/(Q-q): fi:
          od:od:  
       end;
        debug(Fract);
        Fract(5, 13);
        Fract(77, 200);

pellsolve := proc (D::posint) local P, Q, a, A, B, i; if type(sqrt(D), integer) then error "D must be a nonsquare integer" end if; P := 0; Q := 1; a := floor(sqrt(D)); A := 1, a; B := 0, 1; for i do P := a*Q-P; Q := (D-P^2)/Q; a := floor((P+sqrt(D))/Q); A := A[2], a*A[2]+A[1]; B := B[2], a*B[2]+B[1]; if Q = 1 and `mod`(i, 2) = 0 then break end if end do; return A[1], B[1] end proc;
pellsolve(107);
                            962, 93
isolve(x^2-107*y = 13);
       /                             2             \   
      { x = 21 - 107 _Z1, y = 107 _Z1  - 42 _Z1 + 4 }, 
       \                                           /   

         /                             2               \ 
        { x = 86 - 107 _Z1, y = 107 _Z1  - 172 _Z1 + 69 }
         \                                             / 
genpellsolve := proc (D::posint, N::integer) local t, u, L1, L2, sols, x, y; if type(sqrt(D), integer) then error "D must be a nonsquare integer" end if; t, u := pellsolve(D); if 0 < N then L1 := 0; L2 := floor(sqrt((1/2)*N*(t-1)/D)) elif N < 0 then L1 := ceil(sqrt(-N/D)); L2 := floor(sqrt(-(1/2)*N*(t+1)/D)) else return {[0, 0]} end if; sols := {}; for y from L1 to L2 do x := sqrt(N+D*y^2); if type(x, integer) then sols := sols union {[x, y]}; if `mod`(x^2+D*y^2, N) <> 0 or `mod`(2*x*y, N) <> 0 then sols := sols union {[-x, y]} end if end if end do; DEBUG(); return sols end proc;

Debug inefficace