APL(NARS), 35 chars

{⎕FPC←4×⍵+4⋄⍵⍕+/3√1r2+1 ¯1×√23r108}

If one solve the equation

3X^3-X-1=0

one obtain the solution

X=(0.5+m)^(1/3) + (0.5-m)^(1/3) where m=sqrt(23/108).

inplementing that in APL would be the above function... Test:

  ⍪{⎕FPC←4×⍵+4⋄⍵⍕+/3√1r2+1 ¯1×√23r108}¨⍳10
1.3          
1.32         
1.325        
1.3247       
1.32472      
1.324718     
1.3247180    
1.32471796   
1.324717957  
1.3247179572

with last digit approssimate to the remain of the number in memory. This would be the first 300 decimal digits of the Plastic number with last digit approssimate for the remain of number in memory.

  {⎕FPC←4×⍵+4⋄⍵⍕+/3√1r2+1 ¯1×√23r108}300
1.324717957244746025960908854478097340734404056901733364534015050302827
  851245547594054699347981787280329910920994742207425108902639045897
  795594314757096723471754166839038867418751736931584253549908246622
  354533727350458987990956815062774550980248621301216989415752457454
  862507562652461036893890483993227

Bespoke, 197 bytes

cube X=plus X with I;whats X
Richard Padovan,Axel Thue&GH Hardy knew this as ratio
solution came by square-dividing
solution is also morphic,and so^s number X=phi
so this"plastic number"is called P

Calculates the Padovan sequence \$P\$, and outputs \$\frac{P(n+3)}{P(n+2)}\$ as a fraction. (It loops forever if given an input of 0, but that's okay.)

AWK, 34 bytes

{for(;$1^3>$1+1;)$1=($1+1)^(1/3)}1

Step by step:

{                  # $1 is the input (a positive integer)
for(;$1^3>$1+1;)   # as long as $1 cubed is bigger than the input plus one
  $1=($1+1)^(1/3)  # assigns $1 this new value: cube root of $1 plus one
}
1                  # prints $1

Try it online!

Rust, 66 bytes

|n|{let mut i=(1,1,1);for _ in 0..n{i=(i.1,i.2,i.0+i.1)}(i.1,i.2)}

Try it online!

Uses the padovan sequence to approximate the plastic number. Returns a tupled (num,denom). A readable version that works with arbitrary precision integers can be found here.

APL (Dyalog Extended), 19 bytes

{{(⍺+⍵*÷3)*÷3}/⍵/1}

Try it online!

Uses the cube root approximation method, which very, very slowly creeps towards the value of the plastic number.

Explanation

{{(⍺+⍵*÷3)*÷3}/⍵/1} ⍵ → input
               ⍵/1} replicate 1 ⍵ times
 {           }/     reduce by inner fn:
          *÷3       cube root of
     ⍵*÷3           cube root of right arg
   ⍺+               plus left arg

Mathematica, 27 bytes

Nest[(1+#)^(1/3)&,1,#]~N~#&

Uses a truncated approximation of the nested cubic radical form ³√(1+³√(1+³√(1+...))). While the output will always have x-1 decimal places, the result is actually less accurate than that, because the expression converges more slowly than one digit per iteration (x is also used as the number of nested radicals that are computed). For example x = 100 gives

_________________________________________________________________________
1.324717957244746025960908854478097340734404056901733364534015050302827850993693624204577670741656151

where the overlined part is correct.

Ruby, 35 bytes

->x{n=1;n+=0.1**x while n**3-n<1;n}

Try it online!

NewStack, 14 bytes

¹Fᵢ{E2x³⁺÷3x²⁻

Break down:

¹                Add arbitrary number 1 to the stack.
 Fᵢ{             Define for loop with a user's input amount of itterations.
    E            Define new edit for element 0 (element 0 being the 1 added. earlier).
     2x³⁺÷3x²⁻   update x to equal (2x^3+1)/(3x^2-1). (x = element 0).

How it works:

The formula (2x³+1)/(3x²-1) comes from the simplification of Newton's method for the equasion x³=x+1. You can find it here. Repeating this process an infinite amoune of times converges to the plastic number. It's rate of convergence is rather quick at around 2.6 decimals per iteration.

INPUT ITERATION >> VALUE
0 >> 1
1 >> 1.5
2 >> 1.3478260869565217
3 >> 1.325200398950907
4 >> 1.3247181739990537
5 >> 1.3247179572447898
6 >> 1.324717957244746    <- 16 decimal precision in 6 iterations!
...
100 >> 1.324717957244746

Padovan sequence alternative, 27 25 17 bytes

¹Fᵢ{[ƨ2+ƨ3]ℲƤƨ/ƨ2

Break down:

¹                  Append first element of Padovan sequence.
 Fᵢ{       Ⅎ       Define for loop of user's input amount of iterations.
    [ƨ2+ƨ3]        Append sum second and third to last elements.
            Ƥƨ/ƨ2  Print ratio of last two elements.

-2 bytes by choosing better print strategy

-8 bytes by choosing better way to index stack

How it works:

As the Padovan sequence continues, the ratio of the last two elements converge to the plastic number.

INPUT ITERATION >> VALUE
0 >> 1
1 >> 2
...
10 >> 1.3157894736842106
...
89 >> 1.324717957244746    <- 16 decimal precision in 89 iterations
...
100> > 1.324717957244746

Charcoal, 28 bytes

ＡＩθθＡθνＷ‹∕∕Ｘν³θθ⁺νθＡ⁺ν¹νＩ∕νθ

Try it online! Link to verbose mode. Also I apparently messed up Divide and IntDivide :|
Uses the same method as the Python and JavaScript answers.

Axiom, 96 bytes

h(n:NNI):Float==(n>1.E5=>-1;n:=n+1;j:=digits(n::PI);r:=solve(x^3-x=1,10.^-n);digits(j);rhs(r.1))

results

(31) -> [h(i) for i in 0..10]
   (31)
   [1.0, 1.3, 1.33, 1.325, 1.3247, 1.32472, 1.324718, 1.324718, 1.32471796,
    1.324717957, 1.3247179572]
                                                         Type: List Float

how you can see h(2) should be 1.32 and not 1.33 so there is some error in last digits

Then there would be this one of 110 bytes

g(n:NNI):Float==(n>1.E5=>-1;n:=n+1;j:=digits(n::PI);x:=sqrt(23./108);r:=(.5+x)^(1/3)+(.5-x)^(1/3);digits(j);r)

It use the formula for resolve equation of III grade of type x^3-3*p*x-2*q=0 in the case q^2-p^3>=0 that is m=sqrt(q^2-p^3) and x=(q+m)^(1/3)+(q-m)^(1/3)

In our case r^3-r-1=0 this can be written as r^3-3*(1/3)r-2*(1/2)=0 so p=1/3 q=1/2 m=1/4-1/27=23/108 x=(0.5+m)^(1/3)+(0.5-m)^(1/3)

this one that use Newton iteration with start point r=1

f(n:NNI):Float==(n>1.E5=>-1;n:=n+1;j:=digits(n::PI);e:=10^-n;r:=1.;repeat(v:=(r^3-r-1)/(3*r^2-1);abs(v)<e=>break;r:=r-v);digits(j);r)

it change in the function, digits value for obtain one obj of n+1 digits afther the float point. At end the digits() value is assigned back to preciding value.

M, 15 14 bytes

²×3’
*3Ḥ‘÷Ç
Ç¡

Try it online!

Algorithm

This uses rationals and Newton's method. Specifically, for input x, the first x iterations with starting value x are applied.

We're trying to find a specific root of the polynomial p(t) = t³ - t - 1. Newton's method achieves this by taking a starting value t₀ – sufficiently close to ρ – and recursively defining a sequence by
t_n+1 = t_n - p(t_n) / p'(t_n).

Since p'(t) = 3t² -1, we get
t_n+1 = t_n - (t_n³ - t_n - 1)/(3t_n² - 1) = (3t_n³ - t_n - t_n³ + t_n + 1) / (3t_n² - 1) = (2t_n³ + 1) / (3t_n² - 1).

Note that the initial approximation x gets progressively worse as x increases. While the output for x = 3 is slightly less precise than the output for x = 2, since Newton's method converges quadratically to ρ, this shouldn't be an issue for large values of x.

How it works

Ç¡    Main link. Argument: x

Ç¡    Call the second helper link x times, which initial argument x.


*3Ḥ‘÷Ç  Second helper link. Argument: t

*3      Compute t³.
  Ḥ     Unhalve; yield 2t³.
   ‘    Increment; yield 2t³+1.
     Ç  Call the first helper link with argument t.
    ÷   Divide the left result by the right one.


²×3’    First helper link. Argument: t

²       Compute t².
 ×3     Compute 3t².
   ’    Decrement; yield 3t²-1.

PHP, 86 bytes

for($p=[1,1,1];$i++<$argn;)$p[]=bcadd($p[$i],$p[$i-1]);echo bcdiv($p[$i+1],$p[$i],$i);

PHP Sandbox Online

Creates the Padovan Spiral and print the ratio of the last two numbers.

sed, 67 60 (59+1) bytes

s,^,1/1/1 ,
:;s,(1*/(1*)/(1*).*)1$,\2\3/\1,
t
s,(/1*).*,\1,

Try it online!

+1 for the -E flag (ERE instead of BRE). Input and output are both unary: input 11111 for x=5 e.g. Output is a fraction of two unary numbers: the aforementioned 11111 input yields output 11111/1111 (5/4 in decimal).

Approximates the plastic number as a fraction between to consecutive elements of the Padovan sequence.

C#, 317 bytes

using m=System.Math;a=x=>{if(a==0)return "1/1";var d=a(x-1).Split('/');var b=int.Parse(d[0]);var c=int.Parse(d[1]);return string.Format("{0}/{1}",(2*m.Pow(b,3)+m.Pow(c,3)).ToString(new string('#',int.MaxValue.ToString().Length)),(3*m.Pow(b,2)*c-m.Pow(c,3)).ToString(new string('#',int.MaxValue.ToString().Length)));};

It returns the result as a fraction.

Explanation

It uses the Newton's method with x iterations for finding the root of the polynomial p^3-p-1=0. The formula is x_n=1-(f(x_(n-1)))/(f'(x_(n-1))), and x_0 is a starting point.

The polynomials derivative is 3p^2-1, and let's say x_(n-1)=b/c. Then, by using the above formula we get, that x_n=(2 b^3+c^3)/(3 b^2 c-c^3). Let's also say, that we start from 1, this will happen, when x=2, because x>1, and is an integer. Idented, and commented code:

using System;
string PlasticNumber(int x)
{
    if (x == 2) 
        return "1/1";                 

//If x=2, we return our starting value, but we need to return it as a fraction

    var d = PlasticNumber(x - 1).Split('/');
    var b = System.Convert.ToInt32(d[0]);
    var c = int.Parse(d[1]);

//We parse the previous value of the fraction, and put it into two variables

    return string.Format("{0}/{1}", 
        (2 * Math.Pow(b, 3) + Math.Pow(c, 3))
        .ToString(new string('#', int.MaxValue.ToString().Length)),
        (3 * Math.Pow(b, 2) * c - Math.Pow(c, 3))
        .ToString(new string('#', int.MaxValue.ToString().Length)));

//We return the result as a fraction, but it's important not to return it in
  scientific notation, because that will cause issues in the parsing process 

}

MATL (27 28 bytes)

7BG:"t@)y@Q)+h]tG3+)yG2+)/

My first solution (27 bytes)

Try it online!

It's certainly not optimal, I'm still getting used to MATL.

Explanation:

I create a Padovan sequence up to input+3 then find the ratio of the last two numbers.

7B     % Turn 7 into binary to give 1 1 1 
G:"    % For k=1:input do...
t@)    % Existing sequence member k
y@1+)  % Existing sequence member k+1
+h     % Add them together and concatenate to the sequence array
]      % End loop
tG3+)  % Final sequence member
yG2+)  % Second last sequence member
/      % Divide to approximate ρ

Proper fraction output (35 bytes) (28 bytes, @Sanchises):

However, the first solution doesn't fulfill the need for arbitrary precision being the floating point limit of default MATL settings. So rather than adding several bytes to extend this precision, it's simpler to take the proper fraction route and write a fraction of the final two integers in the (N-1)^th and N^th elements of the truncated Padovan sequence.

e.g "114/86"

7BG:"t@)y@1+)+h]tG3+)V'/'YcyG2+)VYc

7BG:"t@tQh)sh]tJ)V47hyJq)Vh&

Courtesy of user @Sanchises. :)

Try it online!

Non-iterative evaluation:

Notably, my shortest code for the 'exact' version is (23 bytes):

1-1h69X^12**108+1I/^6/s

Try it online!

...but doesn't give arbitrary precision. I wonder if anyone can adjust this to fulfill the rules (use the input etc) and still add less than 5 bytes? :P

Octave, 50 bytes

@(n)char(digits(n)*0+vpasolve(sym('r^3-r-1'))(1));

Try it online!

Defines an anonymous function, withn the desired number of digits of output.

This answer abuses that digits returns the current setting for the number of digits in variable precision arithmetic. This means we can just use it in an anonymous function without errors about 'Too many output arguments'.

Other than that, it's really straightforward: vpasolve is short for Variable-Precision Arithmetic Solve, with the precision set by the last call of digits. Since vpa is a Symbolic data type in Octave, which is banned per the spec, we just wrap the whole function in char(...) to get string output. Note that in solve and vpasolve, the f==0 is implied, so r^3==r+1 has been replaced by r^3-r-1 (==0)

05AB1E, 23 bytes

U[X3m¹/¹/¹X+›#X>U]X'/¹J

Try it online!

Direct port of https://codegolf.stackexchange.com/a/126822/59376 by xnor.

TI-BASIC, 21 bytes

:Prompt X //Prompt for input, 3 bytes
:While X  //While X, 3 bytes
:³√(1+Y→Y //Calculate cube root of 1+Y and store to Y, 7 bytes
:DS<(X,0  //Decrement X and skip next command (which doesn't do anything), 5 bytes
:End      //End While loop, 2 bytes
:Y        //Display Y, 1 byte

Uses this recursive formula.

Interestingly, hard-coding the number and rounding it gives the same byte-count:

TI-BASIC, 21 bytes

:Prompt X    //Prompt for input, 3 bytes
:.5√(3       //Store √(3)/2 to Ans, 5 bytes
:Ansֿ¹cosh(3ֿ¹coshֿ¹(3Ans //Store the plastic number to Ans, 9 bytes
:round(Ans,X //Round the plastic number X decimal digits, 4 bytes

Uses this trigonometric formula.

k, 27 bytes

{"/"/$1_|x({y,z,x+y}.)/3#1}

Try it online! This assumes infinite integers (which, alas, is not true). It uses the Padovan sequence.

><>, 38+3 = 41 bytes

11\n;
?!\}2,:01{::::}**-+0(?$-{+{1-:}

Expects the input to be present on the stack at program start, so +3 bytes for the -v flag.

Try it online!

Effectively performs a binary search to narrow in on the output value. Increasing x increases the number of iterations to perform.

Edit: refactored calculation slightly to save 1 byte, previous version:

11\n;
?!\}2,:0{::::}**$-1-0)?$-{+{1-:}

Mathematica, 20 bytes

#^3-#-1&~Root~1~N~#&

Mathematica's builtin Root function gives the solutions to a polynomial equation f[x] == 0.

Explanation

#^3-#-1&~Root~1~N~#&
                   &  (* Function *)
#^3-#-1&              (* A pure-function polynomial, x^3-x-1 *)
        ~Root~1       (* Find the first root *)
               ~N~#   (* approximate to (input) digits *)

Sample I/O

In[1]:= f=#^3-#-1&~Root~1~N~#&;
        f[1]

Out[1]= 1.

In[2]:= f[9]

Out[2]= 1.32471796

In[3]:= f[100]

Out[3]= 1.324717957244746025960908854478097340734404056901733364534015050302827851245547594054699347981787280

Javascript, 36 bytes

f=(x,n=x)=>n**3/x/x<n+x?f(x,++n):n/x

Works the same as the top python answer. No console.log was included because if you run f(x) in console it will be logged automatically.

f=(x,n=x)=>n**3/x/x<n+x?f(x,++n):n/x
console.log(f(300))

Clojure, 46 bytes

#(nth(iterate(fn[i](Math/pow(inc i)(/ 3)))1)%)

Uses the iterated cube-root formula. This is a bit more interesting but longer:

(def f #(apply comp(repeat %(fn[i](Math/pow(inc i)(/ 3))))))

((f 10)1)
1.3247179361449652

Julia 0.5,  44  40 bytes

|(x,r=1)=x<1?r:~-x|big(2r^3+1)//(3r^2-1)

Uses rationals and Newton's method.

Try it online!

Mathematica, 27 bytes

x/.Solve[x^3==x+1>2,x]~N~#&

-1 byte from Martin
-2 bytes from ovs

input

[27]

output

{1.32471795724474602596090885}

Python 2, 49 bytes

n=x=input()
while n**3/x/x<n+x:n+=1
print n,'/',x

Try it online!

The idea is to express the ρ with ρ³=ρ+1 as a fraction n/x whose denominator x is the input accuracy parameter. We take (n/x)³=n/x+1 and clear denominators to get n³=x²(x+n).

Since the LHS increases in n faster than the RHS, we can approximate the equality point n as the smallest with n³≥x²(x+n). The code counts up n until this is the case, starting at x which is smaller.

A small byte save is to divide both sides by x² to write n³/x²≥x+n (negated in the while condition). This is floor division in the code, but the fractional part lost is negligible.

A same-length alternative instead puts x as the numerator:

Python 2, 49 bytes

n=x=input()
while x**3/n/n<n+x:n-=1
print x,'/',n

Try it online!